﻿namespace NumericalData
open System
open System.Threading
open NumericalData.Exceptions
open NumericalData
open ExtensionMethods


    ///// {=========================MatrixGenerators=========================}
type MatrixGenerators< 'a when 'a : ( new: unit -> 'a ) and 'a :> Matrix  > = // and 'a :( static member (@/) : 'a -> double -> unit)
        class
            ///-----------------------------Fields-----------------------------
            ///-----------------------------Constructors-----------------------------
            ///-----------------------------Properties-----------------------------
            ///-----------------------------Operators-----------------------------
            ///-----------------------------Methods-----------------------------
            static member AsyncMatrixGenerate (dim:int) (avg:double) (genMethod:Func<int,double,'a>)
                (asyncSuccess:Action<'a>) (asyncFail:Action<exn>) (asyncCancel:Action<OperationCanceledException>) = 
                let successClearFunc    = fun (ans:'a) ->
                                                asyncSuccess.Invoke(ans)
                let failClearFunc       = fun (exc:exn) ->
                                                asyncFail.Invoke(exc)
                let cancelClearFunc     = fun (exc:OperationCanceledException) ->
                                                asyncCancel.Invoke(exc)    

                let asyncGenerate (dim:int, avg:double) = 
                    Async.FromBeginEnd(dim,avg, genMethod.BeginInvoke,genMethod.EndInvoke)

                let asyncOp = 
                    Async.TryCancelled(
                        (async{
                            let! res = asyncGenerate (dim, avg)
                            return res
                        }),
                        cancelClearFunc)
                let cancellationSource = new CancellationTokenSource()                                                                
                Async.StartWithContinuations(   asyncOp,
                                                successClearFunc,
                                                failClearFunc,
                                                cancelClearFunc,
                                                cancellationSource.Token
                                            )
                cancellationSource

            static member private MatrixGenerator (dim:int) (avg:double) (elemGenerator:int->int->Random->double->double) = 
                let ran = new System.Random()
                let matr = new 'a()
                matr.InnerArray <- 
                    [|
                        for i in 0..(dim - 1) do
                            yield [| 
                                        for j in 0..(dim - 1) do
                                        yield (elemGenerator i j ran avg)
                                  |]
                    |]
                matr

            static member IdentityMatrix (dim:int) = 
                MatrixGenerators<'a>.MatrixGenerator dim 0. ( fun i j ran avg -> 
                    if i = j then
                        1.
                    else
                        0. )

            static member GenerateSimpleMatrix (dim:int) (avg:double) = 
                MatrixGenerators<'a>.MatrixGenerator dim avg ( fun i j ran avg -> (double)(avg*2.*ran.NextDouble() - avg) )

            static member GenerateZeroEnterMatrix (dim:int) (avg:double) = 
                MatrixGenerators<'a>.MatrixGenerator dim avg ( fun i j ran avg -> 
                    let var = ran.NextDouble()
                    if var > 0.8 then
                        0.
                    else
                        (double)(avg*2.*ran.NextDouble() - avg) )

            static member GenerateZeroEnterOnDiagonalMatrix(dim:int) (avg:double) = 
                MatrixGenerators<'a>.MatrixGenerator dim avg ( fun i j ran avg -> 
                    if i = j then
                        let var = ran.NextDouble()
                        if var > 0.8 then
                            0.
                        else
                            (double)(avg*2.*var - avg)
                    else
                        (double)(avg*2.*ran.NextDouble() - avg) )

            static member GenerateDiagonalMatrix (dim:int) (avg:double) = 
                MatrixGenerators<'a>.MatrixGenerator dim avg ( fun i j ran avg -> 
                    if i <> j then
                        0.
                    else
                        (double)(avg*2.*ran.NextDouble() - avg) )

            static member GeneratePositiveDiagonalMatrix (dim:int) (avg:double) = 
                MatrixGenerators<'a>.MatrixGenerator dim avg (fun i j ran avg ->
                    if i <> j then
                        0.
                    else
                        ran.NextDouble()*avg + 0.1 ) 

            static member GenerateLinearFreeMatrix (dim:int) (avg:double) = 
                let matr = MatrixGenerators<'a>.IdentityMatrix dim
                let ran  = new Random()
                for i in 0..(dim - 1) do
                    let row = matr.GetRowVector(i)
                    for j in 0..(2*dim) do 
                        let changeIndex = int (ran.NextDouble()*(float dim))
                        let changeKoef = (avg*ran.NextDouble() - avg/2.)/(float dim)
                        let addition = (matr.[changeIndex] + ( changeKoef*row))
                        matr.[changeIndex] <- addition
                matr

            static member GenerateDominanceDiagonalLinearFreeMatrix (dim:int) (avg:double) = 
                let matr = MatrixGenerators<'a>.IdentityMatrix dim
                let ran  = new Random()
                for i in 0..(dim - 1) do
                    let row = matr.GetRowVector(i)
                    for j in 0..(2*dim) do 
                        let changeIndex = int (ran.NextDouble()*(float dim))
                        let changeKoef = (avg*ran.NextDouble())/(float dim)
                        matr.[changeIndex] <- (matr.[changeIndex] + ( changeKoef*row))
                        matr.[i] <- (matr.[i] + (changeKoef*row))
                matr

            static member GeneratePositiveDefinedMatrixWithEigenValues (dim:int) (avg:double) = 
                let LF = MatrixGenerators<'a>.GenerateLinearFreeMatrix dim avg
                let TransponLF = LF.TransponMatrix()
                let diagonal = MatrixGenerators<'a>.GeneratePositiveDiagonalMatrix dim avg
                let eigenValues = diagonal.ToSquareMatrix().DiagonalVector
                let matr = LF.ToSquareMatrix() * diagonal.ToSquareMatrix() *TransponLF.ToSquareMatrix()
                new EigenvaluesMatrixTuple(matr,eigenValues)

            static member GeneratePositiveDefinedMatrix (dim:int) (avg:double) = 
                let LF = MatrixGenerators<'a>.GenerateLinearFreeMatrix dim avg
                let TransponLF = LF.TransponMatrix()
                let diagonal = MatrixGenerators<'a>.GeneratePositiveDiagonalMatrix dim avg
                let matr = LF.ToSquareMatrix() * diagonal.ToSquareMatrix() *TransponLF.ToSquareMatrix()
                matr

//            static member GeneratePositiveDefinedSymmetricMatrix (dim:int) (avg:double) = 
//                let emt = MatrixGenerators<'a>.GeneratePositiveDefinedDiagonalMatrix dim avg
//                let diagonal:SquareMatrix  = emt.Matrix
//                let eigenvalues = emt.EigenValuesVector
//                for i in 0..(dim - 1) do
//                    for j in 0..(i - 1) do
//                        diagonal.[j,i] <- diagonal.[i,j]
//                new EigenvaluesMatrixTuple(diagonal,eigenvalues)
                    
            static member GenerateDiagonalDominateMatrix (dim:int) (avg:double) = 
                MatrixGenerators<'a>.MatrixGenerator dim avg ( fun i j ran avg -> 
                    if i <> j then
                        (double)(avg*2.*ran.NextDouble() - avg) 
                    else
                        (double)(avg*((double)dim)*2.*ran.NextDouble() - avg)  )

            /// Special Diagonal Dominance Matrix Generator
            static member private GenerateMatrixWithDDBasedOnDistributionArray 
                (avg:double) (distributionArray:int[]) = 
                let ran = new System.Random()
                let dim = distributionArray.Length
                let matr = new 'a()
                matr.InnerArray <- 
                    [|
                        for i in 0..(dim - 1)  do
                            let maxIndex = distributionArray.[i]
                            let maxElem = (float dim)*ran.NextDouble()*avg
                            let mutable randomArray = Array.GenerateFilledWithDouble (dim,maxElem)
                            randomArray.[maxIndex] <- maxElem
                            yield randomArray
                    |]
                matr

            static member GenerateMatrixWithImplicitDiagonalDominance (dim:int) (avg:double) = 
                let mixedIndexes = Array.Shuffle [| for i in 0..(dim-1) -> i |] 
                MatrixGenerators<'a>.GenerateMatrixWithDDBasedOnDistributionArray avg mixedIndexes

            static member GenerateMatrixWithDiagonalDominance (dim:int) (avg:double) = 
                let indexes = [| for i in 0..(dim-1) -> i |]
                MatrixGenerators<'a>.GenerateMatrixWithDDBasedOnDistributionArray avg indexes

            static member GeneratePositiveDefinedMatrixWithDiagonalDominance (dim:int) (avg:double) = 
                let LF = MatrixGenerators<'a>.GenerateDominanceDiagonalLinearFreeMatrix dim avg
//                LF.Print()
                let TransponLF = LF.TransponMatrix()
                let diagonal = MatrixGenerators<'a>.GeneratePositiveDiagonalMatrix dim avg
                let matr = LF.ToSquareMatrix() * diagonal.ToSquareMatrix() *TransponLF.ToSquareMatrix()
                matr

            ///-----------------------------InterfaceImplementation------------------
        end
    and 
    ///// {=========================EigenvalueMatrixTuple=========================}
    EigenvaluesMatrixTuple = 
        class
            ///-----------------------------Fields-----------------------------
            val m_eigenvalues:Vector
            val m_matr:SquareMatrix

            ///-----------------------------Constructors-----------------------------
            internal new(matr:SquareMatrix,evs:Vector) = 
                {
                    m_matr = matr
                    m_eigenvalues = evs
                }

            ///-----------------------------Properties-----------------------------
            member x.EigenValuesVector
                with get() = x.m_eigenvalues

            member x.Matrix
                with get() = x.m_matr
            ///-----------------------------Operators-----------------------------
            ///-----------------------------Methods-----------------------------
            ///-----------------------------InterfaceImplementation------------------
        end